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Abstract 

We investigate the dynamics of an effectively one-dimensional Bose-Einstein condensate (BEC) 
with scattering length a subjected to a spatially periodic modulation, a = a(x) = a(x + L). 
This "collisionally inhomogeneous" BEC is described by a Gross-Pitaevskii (GP) equation whose 
nonlinearity coefficient is a periodic function of x. We transform this equation into a GP equation 
with constant coefficient a and an additional effective potential and study a class of extended 
wave solutions of the transformed equation. For weak underlying inhomogeneity, the effective 
potential takes a form resembling a superlattice, and the amplitude dynamics of the solutions of 
the constant-coefficient GP equation obey a nonlinear generalization of the Ince equation. In the 
small-amplitude limit, we use averaging to construct analytical solutions for modulated amplitude 
waves (MAWs), whose stability we subsequently examine using both numerical simulations of the 
original GP equation and fixed-point computations with the MAWs as numerically exact solutions. 
We show that "on-site" solutions, whose maxima correspond to maxima of a(x), are significantly 
more stable than their "off-site" counterparts. 



2 



PACS: 05.45.-a, 03.75.Lm, 05.30. Jp, 05.45.Ac 

Keywords: Bose-Einstein condensates, periodic potentials, multiple-scale perturbation 
theory, Hamiltonian systems 



I. INTRODUCTION 

Among the most fundamental models studied widely in applied mathematics and em- 
ployed for the description of an extremely large variety of physical phenomena are nonlinear 
Schrodinger (NLS) equations [1]. Applications of NLS equations have become even more 
prominent in recent years due to enormous theoretical and experimental progress that has 
taken place in studies of Bose-Einstein condensates (BECs) 2] and nonlinear optics 
Through the formal similarity between coherent matter waves and electromagnetic waves, 
research on these topics has become closely connected, and progress in one area frequently 
also benefits the other. The cross-fertilization between these two fields is extremely impor- 
tant not only from a theoretical perspective but also for applications. For example, coherent 
matter waves can be manipulated using devices such as atom chips |4J whose design was 
initially suggested by previously-developed optical counterparts. 

In the mean-field approximation, and at sufficiently low temperatures, the dynamics of 
matter waves are accurately modeled by a cubic NLS equation incorporating an external 
potential. In this context, it is called the Gross-Pitaevskii (GP) equation [2|, 

h 2 



-V 2 + s|vI/| 2 + V(f) 
2m 



(1) 



where m is the atomic mass, \I/(r, t) is the macroscopic wave function of the condensate, 
normalized to the number of atoms N (so that J \^\ 2 dr = N), V(r) is the external potential, 
and the effective interaction constant is g = (47rh 2 a/m)[l + 0(( 2 )], where a is the s-wave 
scattering length, and ( 2 = | \I/ 1 2 1 « | 3 measures the density of the atomic gas 

The properties of BECs - including their shape, collective nonlinear excitations (such as 
solitons and vortices), and fluctuations above the mean-field level - are determined by the 
sign and magnitude of the scattering length a. Accordingly, one of the key tools used in 
current studies of BECs relies on adjusting a (and hence the nonlinearity coefficient g defined 
above). A well-known way to achieve this goal is to tune an external magnetic field in the 
vicinity of a Feshbach resonance . Alternatively, one can use a Feshbach resonance 
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induced by an optical [8( or dc electric [9j field. In low-dimensional settings, one can also 



10 



Hi- 



tune the effective nonlinearity by changing the BEC's transversal confinement 

Adjusting the scattering length globally (i.e., modifying a in a spatially uniform manner) 
has been crucial to many experimental achievements, including the formation of molecular 
condensates [12] and probing the so-called BEC-BCS crossover [13J. Additionally, recent 
theoretical studies have predicted that a spatially uniform but time-periodic modulation 



of the scattering length, with a (t) periodica 



attractive condensates in two 



141 ] and three 



ly changing its sign, can be used to stabilize 



151 ] dimensions and thus help to create robust 



matter- wave solitons 161 ] . [In the three-dimensional case, the Feshbach resonance-based 
technique should be combined with a quasi-one-dimensional (quasi-lD) periodic potential, 
known as an optical lattice (OL).] Other nonlinear waves besides solitons have also been 
studied in this context. For example, it was predicted that spatially periodic or quasiperiodic 
2D patterns resembling the classical Faraday ripples in hydrodynamics can be excited in a 
BEC via a nonlinear parametric resonance induced by the spatially-uniform and temporally- 
periodic modulation of a (but, in general, without changing the sign of a) 171 ] . 

More recently, the possibility of varying the scattering length locally (i.e., spatially) has 



also been proposed 
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2l| . Such spatial dependence of the scattering length, which 



can be implemented utilizin g a s patially inhomogeneous external magnetic field in the vicin- 
ity of a Feshbach resonance renders the collisional dynamics inhomogeneous across 
the BEC. Condensates with a spatially inhomogeneous nonlinearity have recently attracted 
considerable attention, as they are relevant to many significant applications, including adi- 
abatic compression of matter waves [isl, Q], Bloch oscillations of matter- wave solitons 20], 



atom lasers 
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231 ]. enhancement of transmittivity of matter waves through barriers 
251 ]. and the dynamics of matter waves in the presence of periodic or random spatial 



variations of the scattering length 



261 ] . Additionally, while these works chiefly concentrated 
on quasi-lD condensates, studies in a quasi-two-dimensional (quasi-2D) setting were also re- 
cently reported [27| (see also Ref. 28] for a qualitatively similar model in a different physical 
context). 

In the present work, we consider an effectively one- dimensional BEC whose scattering 
length is subjected to a periodic variation: a(x) = a(x + L) for some period L. We consider 
the case of no zero crossings, so that a(x) always has the same si gn. Such a nonlinear lattice 
can be realized experimentally (see, e.g., the diagram in Ref. 23]), and it offers various 
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possibilities for the study of matter- wave solitons (as discussed in Refs. [26|, |27|, |29|). Here 
we focus on the dynamics of spatially extended states (rather than solitons), which have not 
yet been considered in earlier works in the context of nonlinear lattices. Extended states 
can be created experimentally, as has been done, for example, in the setting of a BEC in an 
optical superlattice [31] . 

The analysis of the problem under consideration proceeds as follows. First, we apply a 
transformation to the quasi- ID GP equation with the nonlinearity coefficient periodically 
modulated in space in order to derive an effective GP equation with a = const. With 
this transformation, the original inhomogeneity of a(x) is mapped into an effective linear 
potential taking the form of a modified superlattice. We consider spatially extended solu- 
tions of the transformed GP equation in the form of modulated amplitude waves (MAWs), 
whose slow-amplitude dynamics we derive using an averaging method. These analytical 
considerations, presented in Section II, are followed in Section III by an investigation of 
the dynamics and stability of the MAWs using both direct numerical simulations of the 
original GP equation and fixed-point computations with the MAWs as numerically exact 
solutions. We thereby show that "on-site" solutions, whose maxima correspond to maxima 
of the nonlinear lattice, tend to be more stable than "off-site" solutions. We summarize our 
findings and present conclusions in Section IV. We then discuss some technical details of the 
averaging procedure in an appendix (Section V). 



II. THE PERTURBED GROSS-PITAEVSKII EQUATION 

The 3D GP equation ([1]) can be reduced to an effectively ID form provided the transverse 
dimensions of the BEC are on the order of its healing length and its longitudinal dimension 



is much larger than the transverse ones 



11 



321 ] . In this regime, the ID approximation 



is derived by averaging Eq. (DO) in the transverse plane. One can readily show (see, e.g., 
Ref. j20]) that the resulting ID GP equation takes the dimensionless form 

iut = ~Uxx + g{x)\u\\ + V(x)u, (2) 

where we recall that nonlinearity coefficient g(x) varies in space. In Eq. (J2J), u is the mean- 
field wave function (with density \u\ 2 measured in units of the peak ID density n ), x 



and t are normalized, respectively, to the healing length £ = h/ Wno\gi\m and £/c (where 
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c = \/no\gi\/m is the Bogoliubov speed of sound), and energy is measured in units of the 
chemical potential \i = giUQ. In the above expressions, g\ = 2fkJ±ao, where u± denotes the 
confining frequency in the transverse direction, and a is a characteristic (constant) value of 
the scattering length relatively close to the Feshbach resonance. Finally, V(x) is the rescaled 
external trapping potential, and the x-dependent nonlinearity is given by g(x) = a(x)/ao, 
where a(x) is the spatially varying scattering length. 

We will examine periodic modulations of the scattering length by assuming 

g{x) =#o + ^o sin 2 (kx) , (3) 

where Vq and k are, respectively, the amplitude and wavenumber of the modulation. In 
experiments, this modulation mode can be induced by a periodically patterned configuration 
of the external (magnetic, optical, or electric) field that controls the Feshbach resonance. 
We chiefly focus on the case of repulsive BECs, with g > and Vq > 0. For attractive 
BECs, for which g < 0, one can also transform the nonlinear partial differential equation 
([2]) into an effective GP equation with a constant nonlinearity coefficient. 



We now introduce the transformation, v = g(x)u, which casts Eq. ([2]) in the form 

iv t = -~v xx + \v\ 2 v + V(x)v + V e s(x)v , 

V cS (x) = - 2 J- — + J^' (4) 

where f(x) = y/ g{x), and /' = df/dx. Obviously, this transformation applies only in the 
case when g(x) does not cross zero. 

If the scattering-length modulation is weak, i.e., Vq go in Eq. ([31), then V^flf can be 
approximated as a superlattice potential (because it contains a second harmonic in addition 
to the fundamental one) plus a first- derivative operator term: 



-V , , 3k 2 V 2 k 2 V , n , 3k 2 V 2 , a . 
V ea {x) = — + — — - cos(2ks) + cos(4ra) + 



16^o 2 9o 1%. 



Y^H. sin(2/tx) 
2^o 



In this case, we define a small parameter, e = tzVo/go, evincing the fact that the cos(4:Kx) 
harmonic in Eq. is a small correction to the fundamental one, cos(2kx). Thus, to order 
O(e), the effective potential (jSJ) is a modified lattice rather than a superlattice. 

It is also worth mentionioning that in the case of a temporal modulation of the scattering 
length, i.e., g = g(t) in Eq. (|2J), the same transformation defined above can be used (i.e., 



v = \J g{t)u) provided g(t) never vanishes. As is well known, the time-dependent coefficient 
in front of the nonlinear term is then translated into a linear dissipative term on the right- 
hand side of the equation for v(x, t) [30(]. This term has the time-dependent coefficient 
—i(2g)~ 1 (dg/dt)v. While this latter setting may be of interest in its own right (perhaps 
especially when the dissipation coefficient is constant so that the time dependence can be 
completely factored out - i.e., for g(t) of an exponential form), we do not pursue it further 
here. We remark that if g(t) is a periodic function, the transformed equation will generate 
patterns with a periodically oscillating amplitude. Finally, one can complicate the situation 
still further by considering simultaneous spatial and temporal modulations [i.e., g = g(x,t)}, 
though such a configuration would be very difficult to controllably implement in experiments. 



III. MODULATED AMPLITUDE WAVES 



We consider solutions to Eqs. (14115!) in the form of coherent structures given by the ansatz 

v(x, t) = R(x) exp [i [9(x) — fit}) . (6) 
Such nonlinear waves, called modulated amplitude waves (MAWs), have been studied previ- 



33 



34 



35 



361 ]. Because 



ously in BEC models with optical-lattice and superlattice potentials 
of the spatial periodicity of the potential in Eq. (j4]), MAWs in the linear limit yield the par- 
ticular case of Bloch waves (of the transformed GP equation). Generic Bloch wave functions 
are quasiperiodic in x; periodic ones lie at edges of Bloch bands. 

Inserting Eq. ([6]) into Eq. (jlj) and equating the real and imaginary components of the 
resulting equation yields an ordinary differential equation governing the spatial amplitude 
dynamics. For standing waves (for which 9{x) = 0), this equation is 

3 



d 2 R dR , . 
~TT - £-r~ sm(2KX) + 
dx z dx 



2/i — ek cos(2kx) e cos(Akx) — 2V(x) 



R - 2R 6 = , (7) 



where 2fi = 2fi + (3/8)e 2 . Equation ([7]), known as a nonlinear generalized Ince equation 



371 ] , is reminiscent of the nonlinear Mathieu equation, but with the parametric force acting 



on both R and R' . In the linear limit, one can study Bloch waves by applying Floquet 
theory to Eq. ([7]). In particular, one can employ the method of harmonic balance j^S], in 
which one inserts a Fourier series expansion into the (linear) generalized Ince equation and 
studies the resulting infinite set of coupled linear algebraic equations satisfied by the Fourier 
coefficients. 
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A. Small- Amplitude Solutions 



Let us now consider Eq. (JTj) in the absence of the external trapping potential V(x). As 
shown by Eqs. (jlj) and (j5J), the effective superlattice potential V e ^(x) is already a confining 
potential for the condensate. Seeking small-amplitude solutions, we employ the scaling 
R = ye/2s, which transforms Eq. (JTj) into 

3 



<i 2 s (is 
— 7 - e sm{2Kx) — + 
csar dx 



5 2 — £K cos(2ko;) — -e 2 cos(4kx) 
8 



s - es* = . 



(8) 



where 5 2 = 2fi (implying that p, is positive). The alternative scaling R = ewj\pl would 
produce a nonlinearity of size 0(e 2 ) rather than 0(e) in Eq. (jHJ) and would lead one even 
deeper into the small- amplitude regime. 
Equation (jHJ) is of the form 



where 



s" + 5 2 s = eFi(s, s', x) + e 2 F 2 (s, s', x) , 



Fi(s, s', x) = ks cos(2kx) + s' sin(2/tx) + s 3 
3s 

F 2 (s, s', x) = — cos(4k;c) . 



(9) 



(10) 



As in Ref. 36j, we consider situations near the 2:1 subharmonic resonance, for which 5 = 
±/c in Eq. (jHJ). Assuming a small "detuning" from the exact resonance, we introduce the 
expansion 

6 = K + e6 1 +e?6 2 + 0(e 3 ), (11) 
which we insert in Eq. (JHJ) to obtain 



s" + K 2 s = eGi(s, s', x) + e 2 G 2 (s, s', x) + O (e 3 ) 



(12) 



where 



Gi (s, s', x) = — 25iks + F\ (s, s', x) , 
G 2 (s, s', x) = (-5f - 25 2 k)s + F 2 (s, s', x) 



When e = 0, Eq. (fl2j) has the solution 



(13) 



S = pCOs(/€X + 0) , 



(14) 



with first derivative s' = —K,psin(K,x + 0). We now use the solution (TH|) and its derivative 
as a starting point to apply the method of variation of parameters to equation f|T2l) . We 
therefore seek a solution to Eq. (I12p of the form 



s(x) = p(x) cos(kx + 0(x)) , s'(x) = — /tp(x) sin(/tx + 0(x)) 



(15) 



where p(x) and 0(x) are slowly varying functions. Differentiating the expression for s(x) 
in Eq. ffT5l) . we enforce the following consistency condition with the expression for s'(x) in 
(USD: 



p cos(kx + 0) = p0 sin(KX + 0) . 
Subsequently inserting the result in Eq. ( fT2l) yields 

p' = — — sin(/tx + 0)Gi (p cos(kx + 0), — «p sin(/tx + 



(16) 



,2 



— — sin(/CE + (p)G 2 (pcos(fi;x + 0), — npsin(nx + 0), x) 
cos(/tx + 0)Gi (pcos(ra + 0), — /tpsin(/€x + 0), x) 



K 

£ 



(17) 



Kp 



np 



cos(/tx + (p)G 2 (pcos(kx + 0), — up sin(/tx + 0), x) 



(18) 



Our objective is to isolate the components of p(x) and 0(x) that vary slowly as compared 
to the fast oscillations of cos(kx) and sin(/tx) and to derive averaged equations governing 
the dynamics of these parts. To commence averaging, we decompose p and into a sum of 
the slowly varying parts and small, rapidly oscillating corrections. That is, 



p = c + ewi{c, <p, x) + e 2 v i(c, <p, x) + 0(e 3 ) , 

= <p + ew 2 (c, <f, x) + e 2 v 2 (c, <p, x) + 0(e 3 ) . 

We then substitute Eqs. ( TT9i) into Eqs. ( TTTl) -( fi~8l) and Taylor-expand to obtain 

1 sin(/tx + <p)Gi(ccos(kx + <p), — Kcsin(/tx + <p), x) 



(19) 



dx K 
dvi 



dx 



+ Li(c,<p,x) 



+ 0(£ 3 ) 



(20) 



2 cos(kx + <p)Gi(ccos(/tx + if), — Kcsin(/tx + <p), x) 



dx KC 
dv 2 



dx 



+ L 2 (c, <p,x) 



+ 0(e 3 ) . 



(21) 



The functions w\ and w 2 in Eqs. (TT9l) - (T2TT) should be chosen so as to eliminate all the rapidly 
oscillating terms at order 0(e). At 0(e 2 ), the functions v\ and v 2 are similarly chosen. 
The terms Li(c,(p,x) and L 2 (c, ip, x) in Eqs. (12"UI) - (|2~T1) depend on w±, w 2 , G\, and G 2 . We 
provide expressions for wi, w 2 and L 1; L 2 in the Appendix. (Note that expressions for V\(x) 
and v 2 (x) are not needed until one tackles the third-order corrections, so we do not display 
them.) The functions G\ and G 2 are expressed in terms of c and ip as follows: 

Gi(c, <p, x) = —25\kccos(kx + <p) + kccos(kx + <p) cos(2kx) 

— Kcsm(nx + tp) sin(2/tx) + c 3 cos 3 (/tx + (p) , 

3 

G2(c, </?, x) = (—51 — 25 2 k) ccos(kx + p) + -ccos(kx + cp) cos(4kx) . 

8 

To first order, the slow evolution equations are 

c' = 0(e 2 ) , v' = z(^-^:)+0(e 2 ), (22) 

generating a circle of nonzero equilibria with amplitude c* = a/85ik/3. The generating 
functions Wi(x) and u^^) that yield these equations are shown in the Appendix. This 
process results in the amplitude function 

R(x) = \ j^ c * cos(kx + ip#) = v(x, t = 0) , (23) 

where the phase shift p* is arbitrary. The corresponding MAW is given by v(x,t) = 
R(x) exp(—ifit), as per Eq. (jBJ). 

Returning to the original GP equation (in the absence of an external potential), we obtain 
the MAW u(x, t) = v(x, t)/ ^/g(x, t). We examine its dynamics and stability with direct 
simulations of Eq. (j2J) with V(x) = 0. Results of the simulations for g Q = 2, V = 0.15, 
K = it/8 (so that e ~ 0.0295), 8\ = 3, 5 2 = 1, and p* = are shown in Fig. [TJ These 
parameter values correspond to a 87 Rb (or a 23 Na) condensate in a trap with transverse 
confining frequency u± = 2n x 1000Hz that contains N ~ 10 3 atoms with a peak ID density 
uq = 10 8 m^ 1 . Here, the spatial unit (i.e., the healing length) is £ = 0.3 /zm (£ = 0.9 /xm), 
and the temporal unit is £/c = 0.14 ms (£/c = 0.3 ms). Figure [1] shows that the MAW 
appears to be stable for long times as a solution of the original GP equation. For larger e, 
however, the wave breaks up into solitary filaments (that appear to be localized around the 
maxima of the nonlinear optical lattice; see the discussion below), as indicated in Fig. [2] (for 
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FIG. 1: (Color online) Dynamical evolution of the sinusoidal MAW of Eq. f|23|) with parameter 
values = for go = 2, k = n/8, Vq = 0.15, <5i = 3, and 82 = 1. The left panel shows the 
space-time contour plot of the square modulus (density) \u\ 2 of the solution, and the right panels 
show four snapshots (spatial profiles) of the spatio-temporal evolution. The dynamics illustrate an 
apparent stability of the solution up through at least t = 2000. 
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FIG. 2: (Color online) Same as Fig. [Q but for Vq = 1.5. As early as t = 200, the MAW solution 
starts breaking up into solitary filaments that appear to be localized around the maxima of the 
nonlinear optical lattice (see the discussion in the text). 

Vq = 1.5). In fact, as we will discuss below, the apparently stable MAWs obtained for small 
Vq are actually weakly unstable, although their lifetimes are very long. The normalized time 
t = 2000 amounts to 280 ms for a Rubidium (and 600 ms for a Sodium) BEC in real time, 
suggesting that these MAWs can nevertheless be observed in experiments. 

One can also incorporate the generating functions into the MAWs and examine the sta- 
bility properties of the resulting refined MAWs using direct numerical simulations. These 
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FIG. 3: (Color online) Same as Fig. [T] (and also with Vq = 0.15), but for the refined prediction of 
the initial MAW solution based on Eq. (|24j) ■ The solution also appears to be dynamically stable. 
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FIG. 4: (Color online) Same as Fig. [21 but with the MAW of Eq. (|24p used as the initial condition 
for the time-evolution of the GP equation. Observe that the refined initial condition leads to a 
delayed initiation of the instability (which now begins at about t ~ 300). 



MAWs are given by 

R(x) 



-(c* + ewi) cos(kx + cp* + ew 2 ) = v(x, t = 0) 



(24) 



instead of Eq. fl23|) . We plot the corresponding space-time plots and spatial profiles at 
various time instants for Vq = 0.15 in Fig. [3] and for Vq = 1.5 in Fig. HI All other parameters 
are the same as before. It is readily observed that the instability in Fig. H] sets in later (at 
t ~ 300) than that in Fig. |5] (at t ~ 200), which was obtained with the same parameter 
values using the sinusoidal MAW of Eq. ( 1231) . 

The MAWs in Eqs. ( |23l) and (124j) were determined to order 0(e), so they ignore the 
effects of the second-harmonic component of the lattice, which arise at 0(e 2 ). To estimate 
the influence of this component, we used the same MAWs as initial conditions in direct 
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simulations of the original GP equations in the presence of an additional external optical 
lattice (OL) potential, V(x) = — (3/16)e 2 cos(4/tx) that exactly cancels out the second- 
harmonic component of the lattice. While differences between the results of the simulations 
with and without the compensating OL are not apparent in comparing space-time plots 
side-by-side, one can observe the slow development of small discrepancies by examining the 
time-evolution of the absolute value of their difference. 
To second order, the equations of slow evolution are 



e 2 



c I „ c 2 



4.V 5l "^ )sin(2 ^ 



+ 0(e 



3\ 



3c 2 \ 2 



p' = e[8 1 --) + e 



1 3c% 51c 4 If' c 2 , 
52 ~ + " 256^ + 4^ [ 5l ~TJ C ° s(2 ^ 



+ 0(e*). 

(25) 

In contrast to the first-order equations ( |22i) . the equilibrium points of Eqs. ( 1251) depend on e, 
which is typical in second-order averaging. In studying the dynamics of slow-flow equations 
produced by such a procedure, one obtains, in general, a complicated bifurcation problem, 



which can be investigated by taking e small but fixed (see, e.g., Ref. 39]). For our purposes, 
we simply note that the only difference between Eqs. (|25|) and the slow equations one would 
obtain at second order without the 0(e 2 ) lattice term in Eq. amounts to constant terms 
at the same order [they result directly from the detuning of 8; see Eq. (lllj) ]. Without 
the second-order lattice term, 82 would not be present, and there would be an additional 
constant term, —8 2 /(2k), that is cancelled out here because of the extra harmonic. To study 
the effects of the second-order lattice systematically, one can examine the dynamics starting 



with the alternative scaling R = ew/y/2 (rather than R = ^s/2s, as adopted above), in 
which case the effect of the nonlinearity also emerges at second order. 

In the absence of resonances, second-order averaging yields slow dynamical equations for 
(c, (p) with c' = O (e 3 ) and ip' < 0, so that a fixed-radius circle is traversed in the clockwise 
direction. From here, one can also consider resonances whose effects appear at 0(e 2 ). 



B. Stability 

In Figs. [2]and|H we observe that the MAW initial conditions, derived for small Vo, break 
down rather quickly for larger Vq. In fact, there is a weak instability even for small Vo, 
although the time it takes for the instability to set in is rather long (beyond t = 2000). 
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Thus, the MAW solutions have a good chance to be observed experimentally for sufficiently 
small scattering- length modulations Vq. (Recall that t = 2000 corresponds to 280 ms for a 
87 Rb BEC and 600 ms for a 23 Na one.) 

To investigate this point further, we performed fixed-point computations using Eqs. (|2"5|) 
and (1241) as starting guesses in order to obtain numerically exact stationary states. We show 
the results for Eq. (|23|) in Fig. [5] using a computational domain containing one period of 
the solution. In this context, we have implemented a variant of our finite-difference method 
for the spatial discretization) that incorporates Floquet theory, as is described in, e.g., 
40| | . Adding an exp(i#) term and its complex conjugate at the appropriate locations of our 
stability matrix, we fill in the bands of continuous spectrum by varying 9 within the interval 
[0, 27r). We used a partition of 200 equidistant points in 9 for the results shown here. We 
find that the configuration is always unstable, even for the small Vo that appeared to be 
stable based on direct numerical simulations. Hence, the apparent stability of Figs. [1] and 
[3] arises only because of the fact that the simulation was performed for a finite time (up 
to t — 2000). As indicated in Fig. [5] (where the spectral planes are plotted in the bottom 
panels), the instability is weak for small Vq but becomes strong for large Vq. 



To find more stable solutions, we note that the transformation v = y/g(x)u suggests that 
for the nonlinear OL, solutions prefer to be centered on-site (in contrast to what occurs for 
linear OLs), so that their maxima coincide with maxima of the lattice rather than minima. 
This is not surprising, as the multiplicative terms (i.e., the ones without the d/dx) of V e g(x) 
act as a regular (super) lattice after the transformation and the minima of g(x) coincide 
with the maxima of these terms in V e s(x). (We note in passing that this difference between 
linear and nonlinear OLs has also recently been highlighted in [4l| from a spectral theory 
perspective for bright soliton solutions of the NLS equation.) We observe additionally in 
Figs. [2] and H] that when the MAWs break up, one obtains solutions that are localized on 
maxima of the nonlinear OL. We can take advantage of this observation by using MAWs 
with different phases tp* as initial wave functions in simulations of the GP equation (121 . 

We show the dynamical evolution of on-site MAWs (for which = 7r/2) in Figs. [6] and [7] 
for Vo = 0.15 and Vq = 1.5, respectively. As observed in Fig. the solution remains stable 
past t = 1000 instead of breaking up far earlier, as was the case for off-site solutions. We 
confirm these observations on stability using fixed-point computations of the same type as 
above (but now for solutions with y?* = vr/2), the results of which are shown in Fig. [HJ We 
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FIG. 5: (Color online) Results of a fixed point iteration to identify the solution given by Eq. (|23l) 
as a numerically exact stationary state. The top left panel shows the L 2 norm (density) of this 
state in our computational domain (containing one period of the solution), and the top right panel 
shows the maximal real part of the most unstable eigenvalue of the configuration. Observe that 
the configuration is always unstable. Hence, the apparent stability of Figs. [T] and [3] arises only 
because of the finite time (up to t = 2000) of the direct GP simulations. The middle left and right 
panels show the solution (solid curve) and nonlinear lattice g{x) (dashed curve) for Vo = 0.15 and 
Vq = 1.5, respectively. The bottom left and right panels show the corresponding spectral planes 
of the solutions and indicate the weak instability of the former and the strong instability of the 
latter. 

observe that the maximum eigenvalue now has a much smaller value and that the ensuing 
instabilities here are much weaker oscillatory ones. An interesting feature that we note in 
passing is the formation of "rings" of such oscillatory instabilities (see Fig. [H]). 



The stability results discussed above can also be considered in light of the recent work 42] 
on the development of instabilities for NLS equations with constant nonlinearity coefficients 
and periodic potentials. The authors of Ref. [42j found (among other results) that small- 
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FIG. 6: (Color online) Dynamical evolution of the on-site solutions with 0* = tt/2 for the initial 
waves of Eq. ([23]) [panels (a) and (b)] and Eq. ([2"3| [panels (c) and (d)]. In both cases, Vq = 0.15. 

amplitude, periodic, standing-wave solutions corresponding to band edges alternate in their 
stability, where upper band edges are modulationally unstable in the attractive case and 



lower band edges are modulationally unstable in the repulsive case. The theory in [42j is 
also applicable to the MAW solutions constructed above, as the resonance relation they 
satisfy guarantees their spatial periodicity. 

To examine the onset of modulational instabilities, we consider small perturbations to 
the MAWs of the form u\ exp(At — ifit). [That is, one linearizes the GP equation ([2]) around 
the MAW solutions.] Such perturbations satisfy the Bogoliubov equations, 



L+p = \q, L_q = \p, 

where the amplitude U\ = p + iq and 

1 d 2 

L- = ~2^2 + + 9(x)\umaw\ 2 ~ H, 
1 d 2 

L+ = ~2~t^ + v ( x ) + 3g(x)\u M Aw\ 2 - H ■ 
Recall that in the absence of a linear optical lattice, V(x) = 0. 
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FIG. 7: (Color online) Same as Fig. [6j except for Vq = 1.5. Observe that the solutions now become 
unstable for times longer than t = 1000, highlighting the much weaker nature of the corresponding 
instability in comparison to its off-site counterpart. The fact that the refined ansatz (|24p now 
seems to become unstable at an earlier time is probably due to the significant change in the profile 
of the (exact) solution for large Vq (see Fig. [8|). 



In Ref. 42], it is shown that a sufficient condition for the onset of the modulational 
instability is for A = to be in the interior of a band of the L + operator. When A = 0, 
Eqs. ( 1271) decouple, so that L + p = L_q = 0. With the sinusoidal solution (1231) . we obtain 
the following Mathieu equations for V{x) = 0: 



ld 2 q 

l d 2 p 
2d^ 



+ 



+ 



4 

3ec? 



EC 

H ) + — - cos(2/tx + 2<p*) 
3ec 2 

H ) H cos(2kx + 2y?*) 



P 











(28) 



where we recall that 2fx = 2// + (3/8)e 2 = k 2 , c* = a/85ik/3, and 5i and <p* are free pa- 
rameters. (Note that the factors of g{x) in ( |27l) cancel out with g(x) factors in |mmaw| 2 
because of the transformation used to construct the MAW solutions.) In the Mathieu equa- 
tion L + p = 0, the parameter combinations are 3ec^/4 = 2e5\n and /i = k 2 /2 — 3e 2 /16. In 
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FIG. 8: (Color online) Same as Fig. \E\ except for on-site MAWs. In this case, the instability of the 
solutions is much weaker and is oscillatory in nature. 



the limit e 



(i.e., Vr 



o 



0), equations ( |28l) become linear harmonic oscillators. 



To examine the applicability of the aforementioned theorem of [42] in the present setting, 
we have computed the eigenvalues of L + . In particular, as a systematic measure for the 
inclusion of the zero eigenvalue in a band of eigenvalues of L+, we consider min { |Al + | }, 
where Xl + denotes the eigenvalues of the operator. We show this diagnostic in Fig. [9] for 
solutions with 0* = 0. One can see that for Vq < 0.5 (roughly), the eigenvalue is, within the 
levels of our numerical accuracy, contained in a band of L + . This, according to the results 
of 42j, implies the existence of modulational instability in this case. On the other hand, 
for larger Vq, this no longer occurs. However, because the criterion only offers a sufficient 
condition, the results are inconclusive in this latter case. Nevertheless, the full numerical 
results of Fig. [5] indicate the presence of instability here as well. For = n/2, we do not 
find the eigenvalue remaining in a band of the L + eigenvalues for increasing Vo- Hence, 
for that case as well, we need to revert to the full numerical results of Fig. [HI 
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One can also examine the sufficient condition of 42| from a theoretical perspective using 
properties of Hill's equations. The eigenvalue problem L + p = \l + P is described by the 
Mathieu equation 



d 2 p 
dx 2 



+ [a + 2/3cos(2x + 20,)]p = O 



where x = KX an d 



1 + 



2A; 



K 



3£ 2 

8k 



— 2 , P 



2e5! 



(29) 



(30) 



From Floquet theory, we construct solutions to (J2"9"|) of the form p(x) = ex P(l(^L+)x)Z{x), 



where Z(x) is a periodic 



(or Floquet) exponent 
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unction with period it and 7 (Al + ) is known as the characteristic 



4A\. 



Expanding Z(x) in a Fourier series and equating each of its coefficients to zero yields an 
infinite-dimensional matrix, whose determinant A (7) (the so-called "Hill determinant") is 
given by 



A( 7 ) = A(0)-- 

sin [fvr^aj 



The Floquet exponents satisfy A (7) = 0, so they are given by 43j, |44j 



7 



2i 

± — sin 

7T 



-1 



A(0) sin" \ — 7T\ . 



1/2- 



(32) 



To determine the spectral bands, one then computes the trace of the fundamental matrix of 
( 1291) . which is given by 4^ 



tr A /(A L+ ) = 2 cos [-iirj(\ L+ )] 



(33) 



The spectral bands of (l29il are defined by the condition |tr^f(Ai + )| < 2. 

Computing A(0) is, in general, nontrivial, but it is permissible to consider a finite- 
dimensional truncation of the (center of the) Hill matrix provided (3 is small. Using the 
example = 0, the five-dimensional truncation of the matrix is given by 
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With the parameter values k = n/8, Vq = 0.15, g = 2, and 5± = 1 (for which (3 = 0.15), 
we obtain the plot of trM(A_L + ) shown in Fig. [TOJ The eigenvalue Xl + = is part of the 
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0.06 




FIG. 9: The minimal-magnitude eigenvalue of L + as a function of Vq for 0* = solutions. The 
value of this diagnostic indicates the presence of modulational instability. 




FIG. 10: (Color online) Trace of the fundamental matrix (vertical axis) versus eigenvalue \l + of 
the L + operator. One obtains spectral bands when the magnitude of the trace is bounded by 2. 
The zero eigenvalue occurs within such a band, indicating that there is a modulational instability. 
The right panel shows a magnification of the left panel. 

band, indicating that a modulational instability occurs. While this semi-analytical approach 
is presented for completeness and is of interest in its own right, its use of a truncated Hill 
matrix and of an approximate analytical solution seems to indicate that it is preferable to 
follow the diagnostic of Fig. presented above. 
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IV. CONCLUSIONS 



In this paper, we have investigated the existence and stability of modulated amplitude 
waves (MAWs) in quasi-one-dimensional Bose-Einstein condensates (BECs) in "nonlinear 
lattices" . In particular, we considered an experimentally feasible situation in which the con- 
densate's s-wave scattering length is modulated periodically in space. Accordingly, we ana- 
lyzed the Gross-Pitaevskii (GP) equation with a spatially periodic nonlinearity coefficient. 
We transformed this GP equation into a new GP equation with a constant nonlinearity coef- 
ficient and (for small modulations of the scattering length) an effective (quasi-) superlattice 
potential. We subsequently studied spatially extended solutions by applying a coherent- 
structure ansatz to the latter equation, which leads to a nonlinear generalized Ince equation 
governing the amplitude's spatial dynamics. In the small-amplitude limit, we used averag- 
ing to construct MAW solutions, whose stability we examined using both direct numerical 
simulations of the original GP equation and fixed-point computations with the MAWs as 
numerically exact solutions. 

While direct simulations suggest stability of the constructed MAW solutions for suffi- 
ciently weak periodic modulations, the fixed-point computations indicate that even those 
are weakly unstable (though they can persist long enough to permit experimental observa- 
tions) for "off-site" solutions, whose maxima do not coincide with maxima of the nonlinear 
lattice. On the other hand, "on-site" MAW solutions are stable for much longer times even 
for large periodic modulations of the scattering length. This may render the latter solutions 
more straightforward to excite and sustain/observe in an experiment in comparison to the 
former (even though both exist for wide regimes of the relevant parameters). The findings 



are consistent with recent theoretical work in 



421 ] on NLS equations with a periodic potential 



and constant nonlinearity coefficient, suggesting that the results reported in that work can 
be extended to models with spatial modulation of the nonlinearity, such as the one studied 
in the present work. 
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V. APPENDIX: FUNCTIONS IN THE SECOND-ORDER AVERAGING 

In this Appendix, we display some of the functions used in the averaging method in the 
main text. The generating functions at first order for the resonant case are 

c I ' (? \ c 

w\(c, <p,x) = — I 5i cos(2/cc + 2p) + — cos(Akx + 2(p) 

2k V Ak J Sk 

c 

+ - cos( Akx + Aip) cos(2/tx) , (34) 

32k 2 Ak 

1 / c 2 \ 1 
w 2 (c, (p,x) = — [ 5i — — sin(2Kx + 2<p) - — sin(4KX + 2ip) 
2k y 2k J ok 

c 2 1 
— sin(4/tx + Aip) — — sin(2Kx) . (35) 

o2k 2 Ak 

The functions L\ and L 2 that appear at second order are given by 

1 c 
Li(c, ip,x)= — — sin(KX + ip) cos(kx + ip)wiDiGi H — sin 2 (/tx + (p)w 2 DiGi 

K K 

+ sin 2 (KX + ip)wiD 2 Gi + csin(KX + cp) cos(kx + (p)w 2 D 2 Gi 

— — cos(kx + <f)w 2 Gi sin(Kx + ip)G 2 (36) 

K K 

L 2 (c,ip,x)= — — cos 2 (/tx + ip)wiDiG 1 H — cos(/tx + ip) su\(kx + (p)w 2 DiGi 

KC K 

+ - cos(kx + ip) sin(Kx + ip)wiD 2 G\ + cos 2 (kx + <p)w 2 D 2 G\ 
c 

+ — sin(«:a; + ip)w 2 Gi H ^ cos(kx + ip)w\Gi cos(kc + (p)G 2 , (37) 

AvC KjC KjC 

where DiGi and D 2 G\ denote, respectively, the derivatives of Gi(s, s',x) with respect to s 
and s'. 

The functions v i and v 2 are obtained by integrating the terms in L\ and L 2 that can be 
averaged out [i.e., the ones that do not appear in Eq. (l2~5l) ]. Note that L 2 has seven terms, 
whereas L\ has only six. The penultimate term in the former is due to the presence of p 
in the denominator of the expression for (f)' and arises at second order in the Taylor series 
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because p = c + ew\. 
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